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Abstract 

We consider two-fluid flow problems in an Arbitrary Lagrangian Eulerian (ALE) frame- 
work. The purpose of this work is twofold. First, we address the problem of the moving 
contact line, namely the line common to the two fluids and the wall. Second, we perform 
a stability analysis in the energy norm for various numerical schemes, taking into account 
the gravity and surface tension effects. 

The problem of the moving contact line is treated with the so-called Generalized Navier 
Boundary Conditions. Owing to these boundary conditions, it is possible to circumvent 
the incompatibility between the classical no-slip boundary condition and the fact that the 
contact line of the interface on the wall is actually moving. 

The energy stability analysis is based in particular on an extension of the Geometry 
Conservation Law (GCL) concept to the case of moving surfaces. This extension is useful 
to study the contribution of the surface tension. 

The theoretical and computational results presented in this paper allow us to pro- 
pose a strategy which offers a good compromise between efficiency, stability and artificial 
diffusion. 



1 Introduction 

A difficult problem in the modelling of two- fluid flows in a bounded domain concerns the 
displacement of the contact line, namely the points which are at the intersection of the 
solid boundary of the domain and the interface separating the two fluids. The difficulty 
comes from the fact that: 

• the interface follows the fluid motion: the normal velocity of a point on the interface 
is the normal velocity of the fluid particle at the same point, 

• the fluid particles near the boundary of the domain tend to have the same velocity 
as the points of the boundary. 

Thus, if the velocity of the points on the boundary of the domain is zero (classical no-slip 
condition for a viscous fluid on a fixed wall), the moving contact line does not move: this 
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is the so-called moving contact line problem. We refer to the review papers [201 [T7] for an 
introduction to this problem. 

We have to keep in mind that the no-slip boundary condition for viscous flows is only an 
approximation of the Navier boundary conditions (which relates the slip velocity relative 
to the wall to the shear rate boundary). The Navier boundary conditions has been assessed 
by molecular dynamics (MD) simulations, but fails in the vicinity of the contact line. One 
might think that continuum mechanics models (at the macroscopic level) are not able to 
represent the movement of the contact line. Nevertheless, hybrid experiments reported 
in |114 fl2| [24] which consists in imposing the slip velocity obtained by MD simulation 
in a continuum model, have shown that continuum models may provide good results. 
Recently, Qian, Wang and Sheng [TBI E] have proposed a boundary conditions, the so- 
called Generalized Navier Boundary Condition (GNBC), which is completly defined at the 
continuum level and which seems to be in good agreement with MD results, including in 
the vicinity of the contact line. We also refer to [18] for a careful study using MD, also 
showing that continuum models can be used to described the dynamics of the contact line. 

In [161 117] . this boundary condition is used with a finite difference scheme and the dis- 
placement of the interface between the two fluids is handled with a phase-field approach. In 
this paper, we show that the GNBC is very natural with a discretization method based on 
a variational formulation, like the finite element method. In addition, the ALE formulation 
we adopt is convenient to study the energy balance and it uses less numerical parameters 
than in a phase-field formulation (for which a free energy for the order parameter needs 
to be introduced, for example). The main drawback of the ALE method compared to 
other methods to follow a moving interface (like volume of fluid methods [131 [9] , level set 
methods \22\ [T9] or phase-field formulation) is that it does not allow very large motion of 
the interface without remeshing, and it does not allow a change of topology of the domains 
occupied by each fluid. On the other hand, it is generally admitted that it is the method 
of choice when a precise computation of the position of the interface is required, with good 
robustness and conservation properties. The ALE formulation has been used in number 
of applications involving two-fluid flows. We refer for example to [7] for an application to 
the modelling of aluminium electrolysis cells, and to Section [6J for an application to flows 
in narrow channels. 

For efficiency and simplicity, many numerical schemes for two-fluid flows decouple the 
fluid resolution and the movement of the interface. The interface displacement is there- 
fore solved explicitly with respect to the computation of the fluid velocity and pressure. 
The variational formulation proposed in this study being well-suited to study the discrete 
energy of the system, we propose to quantify the spurious energy due the gravity and the 
surface tension when the interface displacement is handled explicitly. The Geometric Con- 
servation Law (GCL) is a natural tool to perform this study. To analyse the contribution 
of the surface tension terms, we will have to extend the GCL to surface integral. 

Here is the outline of the paper. In Section [21 the GNBC is introduced. The ALE 
formulation used to implement this boundary condition is detailed in Section [3] We prove 
some lemmata on Geometric Conservation Law in Section [4[ We then study the energy 
conservation properties of the numerical scheme in Section [5] Finally, Section [6J is devoted 
to some numerical experiments which illustrate the theoretical results and investigate some 
alternative schemes. 
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2 The Generalized Navier Boundary Condition 



We are interested in the two-fluid Navier-Stokes equations, posed on a bounded smooth 
domain ^Cf* (with d = 2 or d = 3), and the time interval (0, T): 

— — h div (pu (g) u) — div (77 (Vm + Vit T )) = — Vp + 'yHn^dj] + Pff, 

div(u) = 0, (1) 
dp 

+ div(pix) = 0. 

The equation is posed in the distributional sense. The velocity is denoted by u, the 
density by p, the viscosity by rj and the pressure by p. The vector g denotes the gravity 
acceleration: 

9 = -#e 3 

where (ei, e%, e%) is an orthonormal basis of the physical space. The term ^Hn^Y, is the 
surface tension term that we will describe in detail below. The system is complemented 
by initial conditions (u(t = 0),p(t = 0)). We suppose that p(t = 0) takes two different 
values p± and p 2 - This property is then conserved as time evolves for the function p so 
that each fluid is distinguished from the other by its density. In the following we denote 
by 

Ui(t) = \x e R d , P (t,x) =pi} (2) 

the domain occupied at time t by the fluid i. We suppose that is a smooth domain, 

and we denote by 

E(t) = an 1 (t)ndn 2 (t) (3) 

the interface between the two liquids. The viscosity rj may depend on the fluid, so that 
rj = rj(p) in ([1]). We denote in the following by 

a = 7] (Vw + Vu T ) (4) 

the viscous stress tensor. In the surface tension term ^Hn^Y,, 7 is the surface tension 
coefficient between the two fluids (which is supposed to be constant in the following), n,£ 
is the unit outward vector normal to (see Figure [T]) and H is the mean curvature of 
the interface X positively counted with respect to the normal n^. The distribution <5s is 
defined by: for any smooth function ip 

{fiE,1>) = J 1>dtTE (5) 

where o"s denotes the Lebesgue measure (i.e. the surface measure) on E. 

Let us now describe the boundary conditions. We first suppose the non-penetration 
condition: 

(u — u b ) ■ uqq = on <9f2, (6) 

where uqq denotes the unit outward vector normal to 17 (see Figure [1]) and u b is the 
velocity of the boundary (u — u b is the slip velocity, u b = for a fixed wall). In addition, 
we will suppose for the sake of simplicity that 

u b ■ n an = 0, 

which implies that the the fluid lives in a fixed domain O. 



3 




The fluid being viscous, we have to set another boundary condition on d£l. As ex- 
plained in the introduction, the classical no-slip boundary condition (u = u b on dQ) 
would stick the interface on the wall, which is clearly unacceptable. 

It may be worth recalling that the no-slip boundary condition is in fact an approxima- 
tion of the Navier boundary condition: 

(u - u b ^j ■ r + an dn ■ t = 0, (7) 

where (3 is the slip coefficient. The coefficient (3 is in practice very large, which explains 
that (u — u b ).T = is a good approximation in many cases. 

As explained in the introduction, the Navier boundary condition is not appropriate 
to describe the moving contact line. This is the motivation of the Generalized Navier 
Boundary Conditions (GNBC) introduced in [16] . To state the GNBC, we need to define 
(see Figure H]) the following vectors, defined on the boundary of the interface: tg^, = 
tie x ngn the tangent vector to <9£, m = x and tgn = uqq x tg S . Both sets of 
vectors (tds,n^,m) and (tg^,tQn,nQfi) are positively oriented orthonormal basis. The 
GNBC writes: for any vector r tangent to d£l, 

p(u-u b ^j ■ t + an m ■ t + 7 (m • t dn - cos(0 s )) t m ■ t 5 d z = 0, (8) 

where as above (3 is the slip coefficient, 7 is the surface tension coefficient between the 
two fluids, and 9 S is the static contact angle at the solid surface. The distribution 5qs is 
defined accordingly with ([5]) by: for any smooth function t/j 

{Sss,^) = f i>dla* (9) 

JdY, 

where denotes the Lebesgue measure {i.e. the length measure) on the curve <9£. 
Compared to the classical Navier boundary condition ([7]), we notice the extra term 

(m • t dn - cos(6 s )) , 
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which is called the uncompensated Young stress in [T7]. This quantity measures the dif- 
ference between the dynamic contact angle 9 made by the interface £ and the boundary 
dti (we choose the convention that this angle is measured in the fluid 1, see Figure HD and 
the static contact angle S . In general, the definition of 9 S is a difficult problem both from 
the experimental and the theoretical standpoints. In our framework, 6 S is assumed to be a 
part of the data. The uncompensated Young stress is concentrated along the boundary of 
the interface and vanishes when the dynamic contact angle is equal to the static contact 
angle. 



3 Variational formulation and discretization 

The aim of this section is to derive a variational formulation of the system of equations (pQ) , 
together with the boundary conditions ©-([8]). 

3.1 The weak ALE formulation 

We assume that for any time t > 0, there exists a smooth and bijective mapping At 
from a reference domain Cl (divided into two separate subdomains Oi and Q2 such that 
fl = fl\ U O2) to the current domain such that At(&i) = tii(t) (see Figure [2D- The 
inverse function (with respect to the space variable) of At is denoted A^ 1 . 




Figure 2: The partition of the domain Q and the ALE mapping. 

The velocity of the domain w is defined by: 

d - 

w{t,x) = —A t (x). 



(10) 



For any function ip(t, .) defined on 0, we denote by ip(t, .) the corresponding function 
defined on the reference domain Q by 

i{>(t,x) = iP(t,At(x)). (11) 

For example, the velocity of the domain w on the current frame is defined by 

w(t,x) = w(t,A^ 1 (x)). (12) 

Notice that the functions ip and tp are such that: 



^(t, x) = ^-(t, A t {x)) + w(t, A t {x)) ■ W(t, A t {x)). 



(13) 



The fact that At maps f2j to ^(t) (i = 1 or 2) implies that the velocity of the domain 
satisfies 

w ■ rii = u ■ rii on dQi, (14) 

where i = 1 or i = 2 and n,; denotes the unit outward vector normal to fij. The density /? 
of the fluid is such that: 

p(t,x) = p(At 1 (x)), (15) 

where /S is equal to p\ on f2i and p2 on 

The following functional spaces will be needed, respectively for the velocity u and the 
pressure p: 

V = L 2 (0, T; H* (fi)), M = L 2 (0, T; L 2 (tt)), 

where 
and 



= {it G (^(n))** , u ■ n dn = on dft} , 



L 2 (ft) = [p G L 2 (0), /p = 

We also introduce the test function spaces on the reference domain 

V = W x n (£l), M = L 2 (Cl). 
In the moving frame, the test function spaces are defined by 

V T = {v: [0, T] x n -> M d , »(t, x) = {j(^ _1 (x)), * G V - }, 
M T = {g : [0, T] x fi -> R, q(t, x) = q{^\x)), q G M}. 

Thus, the test functions do not depend on time in the reference frame ft whereas they do 
on the current one: more precisely, let v be in Vr, then for a fixed x G f2, v(t,At(x)) does 
not depend on time while for a fixed xe!l, v(t, x) does. 

We are now in position to state the weak ALE formulation. It is the following coupled 
problem: we look for a function At '■ & — > ^ and (u,p) inV x M such that u(t = 0, .) = Uq 
and: 

• The function At is smooth and maps fij to f2j(i) (i = 1 or 2). The domains Oj(t) 
occupied by each fluid are thus defined by At and the density of the fluid p is defined 
by: 

p(t, x) = p(A t 1 (x)) = pi, for x G 



(16) 



For all (v, q) in Vr x My, 

— /" pu ■ v + f p(u — w) ■ Vu ■ v — ( div(w)pu ■ v 
dt Jn Jn Jn 

+ / - (Vu + Vu T ) : (Vv + Vv T ) - / pdiv(v) 
Jn% Jn 

= —7 / tr(Vsf) daz — /3 / (w — u 6 ) • u 

Jen 

+7 / cos(9 s )t d Q ■ v dlgx + fv, 

/ g div(u) = 0. 
Jn 



(17) 



Of course, At is not uniquely defined by the condition f)16[> : this is the arbitrary feature 
of the ALE scheme. The function At will be precisely defined at the discrete level in 
Section 13.3.31 

Notation: In (|13p . (JTTJ) and in the sequel, the spatial differential operators are taken 
with respect to the Eulerian variable x. We omit to denote this explicitly for conciseness. 
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3.2 Derivation of the weak ALE formulation 

Let us explain how this weak ALE formulation is obtained from the strong formulation (pQ) , 
with the boundary conditions (J!])-©. 

The starting point of this derivation is based on the following Lemma. 

Lemma 1 For any smooth function tp depending on time t and space x, and any smooth 
function (j) such that <b (defined by <b(t,x) = (j>(t,At(x))) is time-independent, we have: 



d f 

— J iP(t,x)<f>(t,x)dx (18) 
dtp 

(bit, x)-—(t, x) + 6(t, x)w(t, x) ■ Vipit, x) + (bit, x)div( wit, x))ib(t, x) dx. 
at 

The first line in (|17p is obtained by multiplying the material derivative in the equation 
on u in ([!]) by the test function v € Vr and integrating over fh 

f dipu) ,. . . f du _ 

/ — - — • v + div( pu ® u) ■ v = / p— ■ v + pu ■ Vu ■ v, 
Jn ot J n dt 

= — I pu ■ v — I pw ■ Vtt • v — div(w)pu ■ v + pu ■ Vu ■ v, 

dt Jn Jn 

where we used successively the equation on p in ([1]) and then (|18p . The weak formulation 
of the terms involving the pressure are classically obtained by integration by parts. It is 
straightforward to obtain the following variational formulation for the term involving the 
viscous stress and the surface tension term: 

/ - (Vu + Vu T ) : (Vv + Vv T ) - / crn dn -v- / ^Hv n E dcr E . (19) 
Jn 2 Jan Js 

We now use the surface divergence formula (see [25] Equation (24) p. 239 or [1], Equa- 
tion (3.8)). For any smooth hypersurface S in M. d (i.e. a submanifold of W 1 with codi- 
mension 1) with a smooth boundary dT, and normal nj^(x) at point x, one has: for any 
smooth function $ : S — ► W 1 , 

- / H3> ■ n s do-Y = / tr(Vs*) da^ - * mdl ds , (20) 

where the surface gradient is defined by: for any smooth vector field X , 

V^X = P^(x)VX, (21) 

where Pe(x) is the orthogonal projector onto the tangent space to £ at point x: 

Py,{x) = Id - nj](x) (g) ns(x). 

Notice that the surface gradient of X only depends on the values of X on the surface S. 
The vector m is the normal vector to dT, in the tangent space of £ pointing outwards of E 
(see Figure [1]). The measure lg s is the Lebesgue measure on dT,. 

Using the surface divergence formula (|20p . the last two terms in (I19j) writes: 



/ aridn-v- / ^Hv ■ n^da^ = - \ (rngn • v + 7 / tr(Vs«) da^ 
Jan Jt, Jan jt, 



7 / v ■ mdlg-£. (22) 

JdT, 
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We now use the GNBC ([5]) to rewrite the first and last terms in the right-hand side of (|22|) : 



- / crn d Q ■ v - 7 / v ■ m dl d T, 
Jan JdT: 

= P (u-u b )-v + -f {m-t m -cos(9 s ))t m -vdl d z-j v-mdl d x, 
Jan JdT JdT 

= (3 (u-u b ) -v - 7 / cos(9 s )t dn • v digs, 
JdQ JdT 

where we have used the fact that v • m = (v • tgo) (m ■ tgn) (since v • ngn = 0). 

With the divergence formula, we eliminate the mean curvature H (which is difficult to 
approximate at the discrete level), and we naturally enforce the GNBC. 

3.3 Discretization 

The discretization is based on a finite element method in space, and an implicit Euler 
time-discretization. The domain f2 = £l l U £l 2 at the beginning of the n-th timestep, 
where Q% is the domain occupied by the fluid i at time i n , plays the role of the reference 
domain = f2i U f^. 

Given the mesh M n = M™ U M% of the domain^ VP = U Q 2 an d the velocity u n 
discretized in a finite element space at time t n , we aim to propagate these two items to 
time t n+ i, using the weak ALE formulation (|17p . 

In addition to (M n ,u n ), let us give ourselves a space discretization of the domain 
velocity w n at time t n . Depending on the scheme (implicit or explicit treatment of the 
displacement of the domain) the velocity w n may depend on u n or on u n+l . We will come 
back to its computation below, in Section [3.3.31 We introduce the application 

A . / (P")i=l,2 ->■ (^r +1 )i=l,2 / 9 oNj 

^ y i — ► x = y + otw [y) 

which might be seen as an approximation of At n+1 ° Af ■ This application defines the 
domain occupied by each fluicH at time t n +±: = A n:n +i(0,f), for i = 1,2. Without 

loss of generality, the time-step St = t n+ i — t n is supposed to be constant. In the sequel, 
our convention is that y denotes a point in (Of)i=i,2 and x a point in {Q™ +l )i=\ y 2- 



3.3.1 Discretization in space 

We consider a finite element discretization of the domain (Of)j=i 2- It is transported by 
the application An^+i to a finite element discretization of the domain (^" +1 )j = i i 2- The 
finite element spaces at time t n for the velocity and the pressure are respectively denoted 
by 

Vh, n c Hjj(fi), M M cL§(n). 

These finite element spaces depend on the time index n, since the mesh is moving. They 
are supposed to satisfy the inf-sup condition (see [U [H] for example) . 

As explained above, we use test functions which follow the deformation of the domain 
given by An,n+i- the test functions at time t n+ i belong to the following spaces: 

V h , n +i = {v(t n+ i, .) : n -> R N ,v(t n+1 ,x) = v(t n ,A~l l+1 (x)),v(t n ,.) G V h>n }, 
1 and therefore the domains occupied by each fluid 

2 A n ,n+i also defines the mesh at time t n +i- for i = 1,2, each node of M" is transported from to 
b y A»,n+ii thus defining the mesh M" +1 of at time t n+1 . 
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M h , n+1 = {q(t n+1 ,.) : n -> R,q(t n+1 ,x) = q(t n ,A n ^ n+1 {x)),q(t n , .) 6 M M }. 

Unless there is a risk of confusion, we omit the index h for the functions belonging to the 
finite element spaces Vh >n or M^n- 

3.3.2 Time discretization and linearization 

We use the following semi- imp licit Euler discretization of fj!7p : for a given u n £ V/^^, 
(n?)i=l,2, ™ n and (f2™ +1 )i=i, 2 , compute (u n+1 ,p n+1 ) G 14, n+ i x M h)ri+1 such that, for all 

(v(t n , .), g(* n , •)) e V M x M M , 

/ pu n+1 ■ v + / p(w n - ™ n ) • Vu n+1 ■ v - [ dw{w n )pu n+1 ■ v 
+ f V«" +1 + (Vu n+1 ) T ) : (Vv + Vv T ) - I p n+1 div(v) 



— / pu n -v-j / tr(V S n+iu)dcr S n+i - /3 I 

Ot Jan Js n + 1 Jd 

+7 / C0s(6 s )t dn -Vdlg^n+l + / • f, 

JdT, n+1 Jn* 



(u n+1 - u b ) 



(24) 



an n+1 



qdiv(u n+l ) = 0. 

0™+! 

Notice that the superscript n (in Q n ) emphasizes that we consider the domain at time t n , 
even if the boundary of the domain is not moving. When we integrate over f2 n+1 , this 
is to indicate that the test functions and the functions p, r\ (whose values are deduced 
from the domains occupied by each fluid) are taken at time t n +\. If a function defined on 
n n appears in an integral over Q n+1 , it means that this function is transported on £l n+1 
by An,n+i- For example, 

/ t divK) u^ 1 -v = P i I div( U " o A£ ) u^ 1 ■ v(t n+1 , .). 

The discretization (|24|) is obtained from the weak ALE formulation (|17|) , In the third 
line of (|24p appear two terms which are required for better stability properties of the 

scheme. The term / — div(tt n ) u n+1 ■ v is standard. It is analogous to the well-known 

modification introduced by Temam of the convective term (see Section III. 5 in [23]) which 
allows to recover at the discrete level the skew-symmetry property of the advection term. 

The second term — / (u n — w n ) ■ «s u n+1 ■ v da where we have used the notation 

2 J^n+l 

Sp = p 2 -pi, (25) 

is in the same vein, but is specific to the context of two-fluid flows. Notice that both these 
terms are strongly consistent: they vanish for the exact solution. They are introduced in 
order to reproduce at the discrete level the energy estimates that can be derived at the 
continuous level (see Section ED. 

The body force is integrated on a domain denoted by 0*. We will consider in the 
sequel different choices: Q* = £l n or £l n+1 or fj n+1 / 2 . By convention, we define the 
quantity J Qn +i/2 99 • v as 

P9-v = - P9-v + - pgv. 

Q«+l/2 l jQn + l I J Q n 
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It is worth noticing that the three choices have the same computational cost as soon as 
the mesh velocity w n is explicitly obtained from the known velocity w n .They nevertheless 
lead to very different behaviors, as shown in Section 16.11 



3.3.3 The complete algorithm 

To complete the presentation of the numerical scheme, it remains to describe how the 
domain velocity w n is computed. The basic requirement is the kinematic condition (|14p . 
which ensures in particular that the nodes of the mesh which are initially on the interface 
remain on the interface. In addition, An,n+l defined from w n by (|23p must be sufficiently 
smooth so that the mesh remains regular enough for finite element computations. 

In the practical problems we are interested in, it is sufficient to adopt the very standard 
method that consists in solving a simple Poisson problem to compute the velocity of the 
mesh. Moreover, we choose the displacement to be along one direction only, so that we 
actually solve a scalar Poisson problem. This choice, which is definitely reasonable in the 
physical situations that we consider, has important favorable consequences on the quality 
of the algorithm. This will be made precise in Section 15.21 In addition, we discretize 
the velocity of the domain w n in space using the same finite element space as for the 
components of u n . 

We may now write the complete algorithm. Suppose that (fif)i=i,2 and (u n ,p n ) are 
known. Then w n , (^™ +1 )j=i,2 and (u n+1 ,p n+1 ) are computed as follows: 

(i) Compute the terms defined on VL n (like — / pu n ■ v dx ] in the system ([2 

V ot J n n J 

(ii) Compute w n = (0,0, w n ). We consider two options: 

— Scheme (Ml) : explicit treatment of the displacement of the interface 



-Aw n 

w n 

dw n 
dn 



0, 



on n?,i = 1,2, 
on £ n , 

on dO,, 



(26) 



— Scheme (M2) : implicit treatment of the displacement of the interface 



-Aw n 

dw n 
dn 



0. 



u n+1 ■ n S n+i 



on f2£" 1 * , i 
on S n+1 , 



1,2, 



E n+1 



on 



dn. 



(27) 



where the components of the normal ns are denoted by (ni,, ni,, riy,). 

(iii) Move the nodes of the mesh according to An.n+i defined by (j23l) , 

(iv) Compute the remaining terms (defined on O n+1 ) in the system (|24D . 

(v) Solve (|24p to determine {u n+1 ,p n+1 ). The resolution is typically performed by a 
GMRES iterative procedure with an ILU preconditioner and (u n ,p n ) as the initial 
guess. 
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In step (ii), the implementation of the Dirichlet boundary condition on w is made easier 
by defining the normals at each node of the discretized surface S n . Such a definition 
is delicate, since S n is only piecewise smooth, and the nodes are typically singular points 
of S n . In practice, following [3], we use approximated normals at each node of the 
interface, by requiring that the Stokes integration by parts formula holds at the discrete 
level. This is one of the ingredient which ensures the exact mass conservation of each fluids 
on the discretized system. We refer to Section 5.1.3.2 in [7] or to [6] for more details. 

Note that scheme (M2) is much more expensive than scheme (Ml) since it requires the 
resolution of a nonlinear system (the velocity u n+1 and the domains being unknown 
when w n is computed). The resulting nonlinear system can be typically solved with a 
Newton algorithm (see for example [21]). For this study, we chose a relaxed fixed-point 
method, solving iteratively systems PH) and (f27|) at each time step. 



4 Surface Geometric Conservation Law 

We establish in this section a few technical results that will be useful to study the discrete 
energy of the system. 

Let <j) be a test function such that <fi (defined by <p(t,x) = (j)(t, At(x))) is time- 
independent. Taking tp = 1 in (|18p . we readily have 



— f 4>{t,x)dx= f (j)(t,x) divw(t,x) dx, (28) 

There is an analogous formula to compute the time derivative of the surface integral 
of (f) (see [2], formula (4.17) p. 355): 

/ <f>(t,.)da Ik = f 0tr(V St H)d(7 St . (29) 



d_ 



The so-called Geometric Conservation Law (GCL) is related to equation (|28p . The 
precise definition of GCL may differ from an author to another (see e.g. [151 [14} \W\ [5]). 
Here we will say that the GCL is satisfied when the equality in (|28|) is preserved after time 
discretization. The following lemma states that the GCL is satisfied as soon as the mesh 
velocity has only one nonvanishing component. 

Lemma 2 (GCL) Suppose that the domain velocity w n has the form (0,0, w n ). Let 4> be 
a function defined on for i = 1 or 2. Then the ALE scheme satisfies the GCL in 

the following sense: 

/ 4>(x)dx- I <f) o A n ,n+\{y) dy 
Jn n+1 Jn n 

= St 4>o A n<n+ i(y)divyW n (y)dy, (30) 
= St I 6(x)div x (w n oA-^ +1 (x)) dx. (31) 

We refer to [6] or to [7] for the proof of this result. Formulas (f30|) and (f3T|) are useful to 
establish the following lemma which will be used to study the effect of the gravity on the 
energy balance. 
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Lemma 3 Under the assumption of Lemma\^ we have: 

I gx 3 Spw n ■ n S n = / gx 3 Spw n n% n , 
Js n jt.™ 

= -TT ( / pgx 3 dx- [ pgx 3 dx)-^-[ 8pg(w n ) 2 n%n, (32) 
ot \Jn™+ 1 Jn™ J 1 Jt,™ 

and 

/ gx 3 Spw n o A~ l n+l -n E n+i = / gx 3 5pw n o A~\ l+l n% n+1 , 

= ~77 ( f pgx 3 dx- f pgx 3 dx)+^[ Spg(w n ) 2 n^„. (33) 

Ot \Jnn + l Jqu J I J^n 

Proof. Formula (|30p applied to the function pgx 3 yields 
/ P9X3 = / pg(V3 + 5tw n {y)) dy + 5t / pg{y 3 + Stw n {y))d y3 w n dy, 

Jn™+ 1 Jn™ Jn™ 

= / pgy3 + St[ / pgw n {y)dy+ / pgy 3 d y . 3 w n dy) + St 2 / pgw n d y . 3 w n dy , 

= / pgys + St f pgd ys {y 3 w n ) dy + ^- f pgd y3 {{w n ) 2 )dy, 
Jn™ Jn™ 1 Jn™ 

= [ pgys-°~t [ Spgy 3 w n n^ n - ^- [ Spg(w n ) 2 n^ n , 
Jn™ Jt,™ * is™ 

from which we deduce (|32p . We recall that 5p is defined by (|25p . For the second formula, 
we use (ED: 



pg%3 = / pg(V3 + 5tw n {y)) dy + 5t / pgx 3 d X3 (w n o A n x n+l {x) J cte, 

= / pg(y 3 + 5tw n (y))dy -St d X3 {pgx 3 ) w n o A~ l n+l {x) dx , 
Jn™ Jn™+ 1 

= pg{y 3 + 5tw n {y))dy -St pgw n o A~ l n+1 {x) dx 
Jn™ Jn™+ X 

-St Spgx 3 n^ n+1 w n o 

= / pgV3 + St( pgw n {y))dy- pgw n o A^^x) dx) 

Jn™ \Jn™ Jn™+ 1 J 

-St Spgx 3 n% n+1 w n o A~ l n+V 
Js™+ 1 

= pgy3-St 2 l pgw n d y3 w n dy) - St Spgx 3 n^ n+1 w n o A~\ +l , 
Jn™ \Jn™ J Jt,™+ 1 

= [ pgV3 -St [ Spgx 3 n% n+1 w n o A~^ n+l - %-( j pgd y3 {w n ) 2 \ , 
Jn™ Je"+ x l \Jn™ J 

= f pgV3 ~St I Spgx 3 nl n+1 w n o A~\ +1 + %-( I 8pg{w n ) 2 ni^\ , 
Jn™ Jy,™+ 1 z \Jt,™ J 

which gives (j33j) . <D 

The above results are useful to study the energy balance at the discrete level in the 
presence of gravity. To study the effect of surface tension on the energy balance, we have 
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to introduce a notion of "Surface GCL", namely we have to check how equality (|29|) is 
preserved by the time scheme. Contrarily to the volume case, we are not able to prove 
that the scheme satisfies a "Surface GCL". Nevertheless, we can establish inequalities 
which will be convenient in the sequel. 

Lemma 4 (Surface GCL) Suppose that the domain velocity w n has the form (0,0,w n ). 
Let 4> be a function defined on T, n+1 . Then, if 8t is sufficiently small so that 



1 + 6t tr (V S n (w n )) > on S n , 
the ALE scheme is such that: 

/ fido-^n+i - / o An, n +l da^n > St (f) o An in+ itr(X7xn(w n )) da^r, 

Likewise, if 5t is sufficiently small so that 



1 -Sttt [Vvn+i (w n oA~\ 



l n,n+l 



> 



on £ 



n+l 



(34) 



(35) 



(36) 



the ALE scheme is such that: 



i (pda^n+i - / 4> ° Ai,n+i daj^n < 5t / 0tr(V S n+i (w n o A n X n+l )) da^n+i . (37) 

Moreover, in both cases, the difference between the two sides of the inequalities A35\) and 
(Wl) is of order 8t 2 in the limit 5t — > 0. 

Proof. We only consider (|35p . the proof of (|37p being similar. By a change of variable 
we have (see Formula (4.9) p. 353]): 



/ (pda^+i = / 4> o A n , n +i |Cof(V^ ni „, + i)n S n| dcrs 

° Ai,n+1 



J 



ns" + 5t | n\-, n d y3 w n — n\\ n d y2 w n 







d0~Y;n, 



where Cof (V^4 nj „+i) denotes the matrix of cofactors of the Jacobian V^4„ in +i. 
Let us consider the difference 



n\, n dy 3 w n — n\ jn d yi w n 
n^n + 5t | n 2 Zn d y3 w n — n\ n d y2 w n 





n\\|2 



\l + 5ttr(V^(w n ))\ 



(38) 



The first term writes: 

1 + 2Ji ((n^n) 2 d ya w n - n\^ n n^ n d yi w n + (n^n) 2 d y3 w n - n|nn|nO y2 u; n ) 
+ 5t 2 (^(n^nd y3 w n - n^nd yi w n ) 2 + (n{i„dy 3 w n - n^„d y2 w n )' 

For the second term, we have: 

tr (Vs" (™ n )) = div(it) n ) - Vit) n n S n • n S n, 

= d y3 w n - n|n f^nf^w/M , 

= n^n (n^nC^w/ 1 — n^nd yi w n ) + n\ n [r^ n .dy 3 w n — nj^nd y2 w n ) 
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The second term in (|38|) thus writes: 

1 + 2St (n|jn (n^nd y3 w n - n\Kid yi w n ) + n|„ (n| n <9 y3 u; n - n\\ n d y2 w n )) 
+ <5i 2 (n^„ (n^n^w™ - rij :n dy 1 w n ) + n| n (n|„d y3 u; n - nlnd^u/ 1 )) 2 . 

Using now Cauchy-Schwarz inequality, we obtain: 

(n|> (ns„9 3/3 w; n - n\\ n d yi w n ) + n| n (n^nd y3 w n - n^dy^ 11 )) 2 

< ((4„) 2 + (n|„) 2 ) ({n^d y3 w n - nl n d yi w n f + (n|„9 ya ^ - n|„^ 2 ^) 2 ) , 



W —Tlx 



t d yi w n ) + {n\\ n d y3 w n — rij : ndy 2 w n ) . 



8w n 



dy 



(39) 



, where C 



This shows that (|38p is nonnegative and bounded from above by CSt 2 
denotes a constant. 

Thus, if St is sufficiently small so that (|34|) is satisfied, we have: 

/ (pdcr^n+i > / 4> o A n , n +i (1 + <fttr (Vgn (io n ))) das™, 

> / (f) O An,n+ld(Txn + 6t (j) o Ai,n+ltr (V S ™ (if™)) cfos™ . 
JE n is™ 

Moreover, if 5t is sufficiently small so that 1 + <5£tr (Ve™ (w n )) > e > on S n , then 
the difference between the two sides of the inequality ([33]) writes: 

0< / cpda^n+i- / (j) o An, n +i dcr^n — St / o Ai,n+itr(VE™Ciu n )) da^, 



/ 

JE' 



+ 1 



n S n + £t I n\ n d y . i w n — n\\ n dy 2 w n 





|1 + (5ttr(Vs"(^ n ))| doE* 



Using now the inequality < |a|— 16| < 2 |b 1 with a = n^+St I n\\ n d y3 w n — n\ n d y2 w n 

\ ' 

and 6 = 1 + St tr(VE™('U' n )), we obtain the following estimate of the difference between 
the two sides of the inequality (|35j) : 



(40) 





/ 0d<j S n+i - / (j) o An, n +i da^n — St / o ^ njn+1 tr(VE"(if n )) da^n, 
Je™+ 1 Je" Je™ 



C 
< — 
e 



dw n 



dy 



L 2 (E") 



L°°(T, r ' 



Yl\ St' 



This concludes the proof. 



Remark 1 By analogy with the GCL, we would say that the scheme satisfies the "Surface 
GCL" if the inequalities in I135\) or \3T\ j were equalities. Looking at the proof of Lemma^ 
(see 139\) ). one can see that the equality in A35\) would be obtained if and only if 

n|„ = 0, 

which is not satisfied in practice. Devising a scheme that allows to recover an equality in 
[35\) or (37\ ) (and keeping equalities in [30]) and [31]) ) is, to our knowledge, an interesting 
open question. 
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5 Energy estimates and time-discretization 



We now investigate the stability in the energy norm of the implicit and explicit numer- 
ical schemes. In Section I5.1|, we first recall the derivation of the energy estimate at the 
continuous level. We then discuss in Section 15.21 to what extent the computation at the 
continuous level can be reproduced on the time-discretized system. In particular, we will 
emphasize the effect of the explicit treatment of the displacement of the free interface on 
the energy balance. We introduce the following notations: 

K ( t ) = o[ p(t,x)\u(t,x)\ 2 (kinetic energy), 
z Jq 

P v (t)= j 7? ^' ^ \Vu(t, x) + (Vu(t, x)) T \ 2 (viscous power), 
Jn 2 

W(t) = / p(t,x)gx3 (potential energy), 
Jn 

\\ u \\dn = \ \ u \ 2 and (u,v)ga = / u v. 
V Jon Jan 

We denote by |£(i)| the measure of the surface namely J* s d<7 E . 
5.1 Energy estimate at the continous level 

Proposition 1 The physical system described in Section satisfies the following energy 
equality: 

dK dW d\E(t)\ „ n . , f 

~T7 I T7 l~ 7 ' ; + Pv + P(u,u- u b ) gu = 7 / cos{6 s )t m -udl 9 x. (41) 
at dt dt Jqy; 

Proof. 

The proof of this proposition is standard, we briefly sketch it for convenience. Multi- 
plying the strong form of the equations by u and integrating, we obtain: 



1 d f . l2 f n , T ,2 , 

- — / + / -|Vu + Vu | + 0{u - u b ,u) dn 

J Q, J Q 

= -7 / tr(V E «) d«T E + / pg u + j cos(9 s )t d n - udlg^. 



(42) 



Here are the main steps to obtain this equation: 



d(pu) f dp. o /" 9m , 

■udx = I — \u\ dx + I pu-—-dx 



dt Jo. dt Jn dt 

Q 



= - f -7T- \u\ 2 dx + -— [ p\u\' 2 dx, 
2j n dt l 1 2dt.A/ 1 1 

= \J Q P U ' V(u 2 )dx + ~ J^p\u\ 2 dx. 

The first term in the right-hand side of the last equality cancels with the term resulting 
from the integration of div(pu <S> u) ■ u. The remaining terms of (|42p comes from the 
integration by parts of the stress terms and the GNBC ([S]). 

We now consider the surface tension term. We have, since w ■ uqq = (see [2[ 
Formula (4.17) p. 355]), 

J tr(V E u)(ia E = J tr(V E u>)d<7 S = -| £ da^. (43) 
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For the last equality, we refer for example to [2], formula (4.17) p. 355. The first equality 
relies on ()20|) and on the fact that u ■ = w ■ «s which implies that u ■ m = w ■ m. 
Concerning the gravity term, we have: 

/ pg u = - pgVx 3 u= div(pu)gx 3 , 
Jn Jn Jn 

-9x3 = --77 / PQX3, (44) 



la at"" dt 

where we have used the fact that Vp = 5p ns (see ([25]) for the definition of dp) . 



5.2 Energy estimates at the discrete level 

In this section, we present the discrete counterpart of the computations previously made 
for the continuous system. For simplicity, the problem is only discretized in time. Nev- 
ertheless, all the following computations could be carried out after space discretization 
with finite elements as soon as the pressure finite element space contains linear functions 
and with a definition of the normals at the interface consistent with the discrete Stokes 
formula (see [HE]). 

We define the discrete counterpart of the notations introduced before: 

K n = - / p\u n \ 2 (kinetic energy), 
2 Jn™ 

P«= [ ^|Vu n + (Vu n ) T | 2 (viscous power), 
Jn™ 2 

W n = / pgx-3 (potential energy), 
Jn™ 

and we denote by |S n | the measure of the surface S n . 

Proposition 2 We first consider the case without body force and without surface tension, 
namely we assume g = and 7 = 0. We suppose that the domain movement is computed 
with the explicit scheme (Ml), i.e. solving equation 126\) . Then, we have the following 
stability result: 

+ PZ +1 + f 3(u n+1 -U b ,U n+1 ) dSl = - / £-\u n+1 oAn,n+l-U n \ 2 . (45) 

dt J Q n 2dt 

The estimate (j45|) is the discrete counterpart of (|4ip . For the proof, we refer the reader to 
[6] where a similar result is presented. It is based in particular on the GCL (130p . Since the 
right-hand side is nonpositive, it shows that the time discretization scheme does not bring 
spurious energy to the system. Such a property is easy to obtain on a fixed mesh, but more 
complicated to ensure on a moving mesh (see [6] or [7] for the details). It is noteworthy 
that such a stability property is obtained with a scheme which only requires to solve a 
linear system at each time step (since scheme (Ml) decouples the mesh movement and the 
fluid resolution). Nevertheless, we have to bear in mind that the assumption g = and 
7 = is not satisfied in the interesting physical configurations. We will show in the sequel 
that when gravity or surface tension are considered, an analogous stability in the energy 
norm can be obtained with an implicit treatment of the domain movement (scheme (M2)) 
but is no more valid with an explicit scheme like (Ml). 
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Proposition 3 We suppose that the domain movement is computed with the explicit 
scheme (Ml), i.e. solving equation Ii26\) and that the gravity is computed on J7 n+1 (namely 
O* = f2 n+1 in (|24|) ). Then we have the following stability result: 

K n+l_ K n W n+2_ W n+l |S"+ 2 | - Ql n+l 

si + si + P " + 7J V [ + ^ - ^ u )d « 

(46) 



with 



and 



--/ <5 W K +1 ) 2 n E „ + i )3 > 0, (47) 

S n + 1 



2 



6^ P = | (V +2 | - - St J tr(V E „ +1 ^ +1 )) > 0. (48) 

Remark 2 In i/ie above Proposition, note that the terms e™ and e™ exp are of order St, 
but are nonnegative: they bring a spurious power to the system. This is a consequence of 
the explicit treatment of the free interface displacement. 

Proof. 

We only focus on the terms related to gravity and surface tension. The other terms 
can be treated as in Proposition [2l following the arguments of [6]. 

Consider first the gravity term. The purpose is to mimic at the discrete level the 
computations done to establish (|44p . Multiplying the equation by u n+1 , we have: 

pg . u n+l = _ f pgVx 3 -U n+ \ 



gx 3 div(pu 



/ gx% 5p u n+1 ■ n E n+i . 



Therefore, noticing that with Scheme (Ml) we have u n+l ■ n S n+i = w n+1 ■ n S n+i , and 
using formula (|32p of Lemma [3] we have: 



/ pgu n+1 = gx 3 5pw n+1 ■ n s „+i, 



S n+1 
1 
St 

'11 

~2 



gxsdpw n Sn+ i, 



pgx 3 dx - I pgx 3 dx 
n n + 2 ,/n n + 1 

Spg(w n+1 ) 2 nl n+1 . (49) 

S n+1 



We note that the last term in (|49p . namely is of order St and it may bring some 

energy in the system since it is nonnegative. Indeed, if the heaviest fluid is below the 
lightest one (which is the case in our practical applications), Sp < and n% n+1 > 0. 
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We next consider the surface tension term and we try to reproduce on the discrete 
system the formula (f4~3]) . Using again u n+1 • n s „+i = w n+l ■ n S n+i, we have with ([20]) . 



-1 



I tr(V S n+itt n+1 ) = -7 / tr(V E n + i^ n+1 ) 



= -1 (|S n+2 | - |S n+1 j) 

+ 2 (|E n+2 | - |S n+1 | - St [ tr(V E n+iK; n+1 )^ . (50) 

To estimate the last term in ((30]) . namely e~j^L, we need the Surface GCL result proved 
in Lemma (U By choosing </> = 1 in (|35p . we see that the last term in (]50p is of order <5t 
and that it is nonnegative in the limit St — ► (see (1401) ). 

Proposition 4 We suppose that the domain movement is computed with the implicit 
scheme (M2), i.e. solving equation A26\) . and that the gravity is computed on Q n+1 (namely 
0* = f2 n+1 in (|24p J. T/ien we /iai>e the following stability result: 

- ~ K + st + K +1 + 7 !S \ t ' £ ' + - «*, W (51) 

+ / " ""I' = ~ e 9 ~ e l^V + 7 / COs(e s )t dn ■ U n+1 dl d x 



with 



and 



e 9=- 6 {l 



-l^ n+1 \-\E n \-Stj^tr(V^w n )^ >0. (53) 



77ie scheme is therefore stable in the energy norm. 
Proof. 

Consider first the gravity term. With the implicit scheme (M2), we have w n ■ n S n- 
u n+i . nsn+1 ( see (|27p ). Thus, using formula (]33]) of Lemma (3] the counterpart of 
reads (compare with (]49D ): 

pg . u n+1 = / gx 3 Sp u n+1 • n s „+i , 

Qn+l Je™+ 1 



1 

'ft 



pgx^dx — / pgx^dx 
n n + L Jn n 



+ S 4 [ Spg(w n ) 2 nl n . (54) 



2 

The last term is exactly e™. Consider now the surface tension term: 

-7 / tr(V s „+i-u n+1 ) = -7 / tr(V E «+i«? re ), 

= _^(|E«+i|_|E»|) 



+ f - - «W / tr(V S n + i«; n ) ) 

St \ / 



(55) 
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To estimate the last term in (f5Uj) . namely — e™ imp , we use again Lemma H] (surface GCL). 
By choosing (j) = 1 in (|37p . we see that it is nonpositive in the limit St — > 0. 



5.3 Discussion 

In summary, we have shown that the numerical scheme with an explicit treatment of the 
free surface displacement is stable in the energy norm in absence of gravity and surface 
tension. But in presence of gravity and/or surface tension, some terms may introduce a 
spurious energy in the system if the treatment of the interface displacement is explicit 
(and if the timestep is not sufficiently small). Moreover, we observe that in both cases 
(gravity and surface tension) the terms of order 5t ( (|49p or (I50j) ) which have the wrong 
sign are integrals over the interface between the two fluids. These theoretical results are 
in agreement with typical observations. We indeed often notice that when a numerical 
instability occurs, it is located on the interface between the two fluids. Moreover, such 
instabilities typically increase with increasing gravity or surface tension. 

We have also shown that an implicit treatment of the free surface displacement yields 
a numerical scheme which is stable in the energy norm. Of course, such a scheme is more 
expensive, and we show by numerical experiments in the next section how to build schemes 
with are cheap but seems to have similar stability properties. 

6 Numerical experiments 

In all the following numerical experiments, the Navier-Stokes equations are discretized 
with the Q2/P1 pair of finite element, with a discontinuous pressure. 

6.1 Energy balance in presence of gravity 

The purpose of the simulations presented in this section is to illustrate on a simple physical 
experiment the theoretical results established in Section We consider two fluids sub- 



Fluid 2 




Fluid 1 



Figure 3: Schematic representation of the two- fluid experiment with gravity. 

jected to a vertical gravity (see Figure [3j). The lowest fluid is the heaviest. The domain 
is fi = (—2,2) x (0,2). The equation of the steady state interface is X2 = 1- At t = 
the interface is perturbated, its equation is X2 = x±/5 + 1. The experiment consists in 
observing how the interface goes back to equilibrium (namely zero velocity, hydrostatic 
pressure, interface X2 = 1). Here are the physical parameters (in reduced units): p\ = 1, 
P2 = 0.91, 771 = 0.01, r]2 = 0.0091, g = 100. In this test case, we neglect the surface tension 
effect (7 = 0) and we assume pure slip on the wall (the GNBC will be investigated in the 
next test case). 
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6.1.1 Explicit treatment of the displacement of the interface 

We first consider the explicit scheme (Ml) with the body force integrated on ft n+1 (namely 
0* = J7 n+1 ). In Figure [U we plot the quantity 



^g,expl 



K n+1 - K n W n+2 - W n+1 
- + — 



+ p: +i + 



St St 
and the dissipation due to the Euler scheme 

_P_ 

» 2St 



fin 25t 



u 



n\2 



c euler 



U 



n+l 



U 



n.|2 



(56) 



(57) 



Two comments are in order. First, we observe that the quantity e™ expl is indeed nonneg- 
ative which is an numerical illustration of the result proved in Proposition [3j Second, we 
observe that the spurious energy provided by the gravity term can be greater than the 
energy dissipated by the Euler scheme. Other simulations (not reported here) confirmed 
that both terms are of order 0(St). 



Euler dissipation 
Spurious power 




Figure 4: Comparison of the spurious power term (|56p and the Euler dissipation (|57p with 
the explicit scheme (Ml). The positiveness of the spurious power are in agreement with 
the results of Proposition [3] and can potentially induce unstabilities. The positiveness of 
the Euler dissipation have of course a stabilizing effect. 

To illustrate the importance of the choice of the domain where the body force is 
integrated, we have represented in Figure [5] the result obtained with Q* = £l n . It is worth 
noticing that the computional cost is the same in both cases (f2* = 0™ or Q n+1 ) since 
the movement of the interface is treated explicitly in scheme (Ml). With St = 0.025 the 
result with 0* = J7 n is similar to the result with ft* = £l n+1 and St = 0.1. But with 
St = 0.1, we observe that the result may be dramatically unstable when f2* = J7 n (we 
indeed observe growing oscillations of the free interface whereas the system is expected to 
go to its stationnary rest state). 



6.1.2 Implicit treatment of the displacement of the interface 

In Figure [H we consider the fully implicit scheme (M2) with 0* = £l n+1 and we plot the 
quantity 



K 



n+l 



St 



K n W n+1 
+ - 



St 



+p. 



n+l 



+ 



fi" 2M 



U 



n|2 



(58) 
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Figure 5: Comparison of the viscous power term P™ when the gravity is integrated on 
Q* = £l n or f2 n+1 (referenced respectively as "explicit gravity" and "implicit gravity"). 
Although both choices have the same computational cost (we use scheme (Ml)), we observe 
that the choice Q* = Q n may induce very severe unstabilities (here when 5t = 0.1). 



and the dissipation due to the Euler scheme. The results are in agreement with Proposi- 
tion the balance is now negative which ensure the stability of the scheme. We neverthe- 
less recall that this interesting feature is obtained to the price of an expensive scheme (since 
a nonlinear problem has to be solved at each time step to achieve an implicit resolution 
of the interface movement). 



6.1.3 Alternative schemes 

In this section we try to see how the theoretical results established so far can lead to simple 
alternatives to the schemes (Ml) and (M2). The purpose is to devise a scheme whose cost 
is the same as (Ml) and which has dissipation properties similar to (M2). 

Comparing equations (j32[) and (|33l) in Lemma El we observe that the last terms have 
opposite signs: in one case it is dissipative, in the other case it brings a spurious energy. It 
is therefore natural to try to combine these two equations in order to decrease the amount 
of spurious energy appearing in explicit schemes. This simple observation leads to propose 
to integrate a part of the gravity on Q n and the other part on Q n . In other words, we 
can take 0* = n n+1 / 2 in (p3[). 

A second natural idea to circumvent the expensive implicit resolution of the free inter- 
face movement is to use an explicit scheme with an extrapolated interface velocity. More 
precisely, we propose to approximate ([59]) with: 

Scheme (M3): displacement of the interface with an extrapolated velocity 

-Aw" = 0, inOf,i = l,2, 

(2u n -u n - 1 oA-\ n )-n Sn 



dn 



on S n , ( 59 ) 



0, on dQ, 



We have tested these two simple ideas in the above experiment (two fluids submitted 
to gravity). The results are reported in Figure [7J First, we set 0* = O n+1 , and we 
compare the scheme (M2) with the scheme (M3). We observe that they have almost the 
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10 

Time 



Figure 6: Comparison of the spurious power term (i58j) and the Euler dissipation (i57|) with 
the (expensive) implicit scheme (M2). Contrarily to the result of Figured] the balance is 
negative, which is in agreement with the results of Proposition |4] and ensure the stability 
of the scheme. 



same dissipation property. Thus, (M3) has (almost) the same stability properties as the 
expensive (M2) scheme and the same computational cost as the "cheap" (Ml) scheme. 
Second, using again scheme (M3), we notice that integrating the gravity on f2* = S7 n+1 / 2 
decreases the artificial dissipation of the scheme. 

Therefore, it seems that a good compromise between efficiency, stability and artificial 
dissipation is to use scheme (M3) with 0* = f2 n+1 / 2 . 

Of course, these properties still have to be assessed in more complex situations. The 
only purpose of this section was to show the potential interest of the energy studies pre- 
sented above to devise new schemes. 

6.2 Moving contact line problems and GNBC 

To illustrate how the GNBC can handle the moving contact line, we present the results of 
two benchmarks proposed in [17] involving a Couette flow for two fluids. The geometry is 
2D and represented in Figure El The domain is periodic along x: the dashed lines, located 
on x = and x = 4L represent the periodic boundaries. The walls are defined by y = 
(bottom) and y = H (top). A velocity Ve x (resp. — Ve x ) is imposed on the top (resp. 
on the bottom) of the domain. For the first test case (the "symmetric" one), we took the 
following values of the parameters from [17] (given in reduced units): H = 13.6, L = 27.2, 
V = 0.25, p! = p 2 = 0.81, 171 = ?? 2 = 1-95, 7 = 5.5, ft = j3 2 = 1.5 and 9 8 = vr/2. The 
parameters for the second test case (the "asymmetric" one) are the same except V = 0.20, 
/?2 = 0.591 and 9 S such that cos 9 S ~ 0.38. In both cases, of course, there is no gravity: 
9 = 0. 

At t = 0, the interfaces separating the two fluids are straight and vertical. After a 
while the interfaces reach a steady state position. Let us emphasize that this behaviour 
is a direct result of the GNBC conditions: the interfaces would of course not converge 
to steady curves if we had imposed u x = V on the top and u x = —V on the bottom 
(no-slip). On the other hand, such a result cannot be obtained with pure slip boundary 
conditions. Figure [9] shows the velocity field at t = 1,10,20,100 in the symmetric case. 
The color represents the magnitude of the velocity. At t = 100 the stationnary state is 
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Figure 7: Spurious power term (|56p for scheme (Ml) (same result as in Figure 0]) and the 
term (|58j) for the implicit scheme (M2), or the explicit scheme (M3) with an extrapolated 
interface (with for all these computations £2* = J7 n+1 ). We observe that when the gravity 
is integrated on £1* = fj n+1 / 2 (with the scheme (M3)), the scheme is still stable (negative 
energy balance) but also less dissipative. 
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Figure 8: Schematic representation of the two-fluid periodic Couette simulation. 
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almost reached: the blue zone surrounding the interfaces shows that they are fixed (which 
means they indeed sleep with respect to the wall), whereas the remaining part of the wall 
is red, which corresponds to a fluid adherence on the wall. More quantitatively, Figure [TU] 
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Figure 9: Velocity field at t=l, 10, 20, 100 in the symmetric case. The color represents 
the velocity magnitude. 

shows the velocity on the top wall at t = 5 (the contact points are still moving) and 
at t = 160 (the contact points are fixed). Figure [TT1 shows the evolution in time of the 
velocity of a point on the contact line (which tends to zero) and of the velocity of a point 
on the wall far from the contact line (which tends to about 0.21, whereas a total adherence 
would correspond to 0.25). The stationary interfaces in the symmetric and asymmetric 
cases are represented on Figure [12j These results are in very good agreement with those 
presented in [17], which were obtained either by a continuum phase-field formulation, or 
by molecular dynamics simulations. 

To illustrate the result established in Proposition^ we performed the above experiment 
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Figure 10: Velocity on the top wall versus x for the symmetric case. At t = 5, we are 
still in the transient phase, the velocity of the points on the contact line is non-zero. 
At t = 160, the points on the contact line are almost fixed (which means they indeed 
move with respect to the wall), whereas the points on the boundary far from the contact 
line move with a velocity of magnitude about 0.21 (0.25 would have meant a complete 
adherence on the wall). 




Figure 11: Velocity of the point CI on the contact line and of the point A on the boundary 
(see Figured]) versus time for the symmetric case. The velocity of the point on the contact 
line tends to zero (pure slip) whereas the velocity of the point A tends to about 0.21. 
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Figure 12: Stationnary interface profiles in the symmetric and asymmetric cases. 



with three values of surface tension 7 = 5.5, 11, 55 and we plot in Figure [T3l the quantity: 
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/ cos{9 s )t an ■ u n+1 dl d x. 



(60) 

According to Proposition [3l this quantity is equal to e"+^ p (defined in ([15]) ). The results 
presented in Figure [13] confirm that the surface tension indeed generates spurious power 
(the balance is positive), which increases with the surface tension coefficient 7. 
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Figure 13: Effect of surface tension on the energy balance (|60f) when the interface dis- 
placement is treated explicitly. Observe that this quantity is indeed positive and increases 
with the surface tension coefficient 7 (at least in the dynamic part of the simulation when 
the interface is still significantly moving). This confirms the results of Proposition [3j 
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7 Conclusion 



We have proposed a formulation for the Generalized Navier Boundary Condition which 
allows to compute the moving contact line at the interface between two fluids. This for- 
mulation is in particular well-suited to an energy stability analysis. We have shown that 
an explicit treatment of the free interface displacement introduces a spurious power in 
presence of gravity and surface tension. The study of the surface tension terms required 
an extension of the Geometric Conservation Law concept to surface integrals. Several 
stability results have been established and illustrated with numerical experiments. Un- 
derstanding the spurious power induced by numerical schemes may help to devise more 
stable algorithms. In particular, we have proposed a simple scheme which offers a good 
compromise between efficiency, stability and artificial diffusion. 

Extension to second order schemes and generalization to more general cases of the 
"Surface Geometric Conservation Law" (Lemma could be investigated in future works. 
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